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We present a numerical method to compute quasiequilibrium configurations of close binary neu- 
tron stars in the pre-coalescing stage. A hydrodynamical treatment is performed under the as- 
sumption that the flow is either rigidly rotating or irrotational. The latter state is technically more 
complicated to treat than the former one (synchronized binary), but is expected to represent fairly 
well the late evolutionary stages of a binary neutron star system. As regards the gravitational field, 
an approximation of general relativity is used, which amounts to solving five of the ten Einstein 
equations (conformally fiat spatial metric). The obtained system of partial differential equations 
is solved by means of a multi-domain spectral method. Two spherical coordinate systems are in- 
troduced, one centered on each star; this results in a precise description of the stellar interiors. 
Thanks to the multi-domain approach, this high precision is extended to the strong field regions. 
The computational domain covers the whole space so that exact boundary conditions are set to 
infinity. Extensive tests of the numerical code are performed, including comparisons with recent an- 
alytical solutions. Finally a constant baryon number sequence (evolutionary sequence) is presented 
in details for a poly tropic equation of state with 7 = 2. 
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I. INTRODUCTION 



Inspiraling neutron star binaries are expected to be among the strongest sources of gravitational radiation that 
could be detected by the interferometric detectors currently under construction (GEO600, LIGO, and VIRGO) or 
in operation (TAMA300). Such binary systems are therefore subject to numerous theoretical studies (see e.g. Q| 
for a review). Among them there are (i) Post-Newtonian (PN) analytical treatments (e.g. p|, jjj) and (ii) fully 
relativistic hydrodynamical treatments, pioneered by the works of Oohara & Nakamura (see e.g. g), Wilson et al. [^J^] 
and recently developed by Shibata [p|-pd|, the Neutron Star Grand Challenge group (l3[R| and Oohara & Nakamura 
Jh|. These last three groups integrate forward in time the evolution equations resulting from the 3+1 formulation 
of general relativity fiq , |l6|| . In parallel of these dynamical calculations, some quasiequilibrium formulation of the 
problem has been developed Jl7|-p0[ and successfully implemented [^l]-|4|. The basic assumption underlying the 
quasiequilibrium calculations is that the timescale of the orbit shrinking is larger than that of the orbital revolution 
in the pre-coalescing state. Consequently the evolution of the binary system can be approximated by a succession 
of exactly circular orbits, hence the name quasiequilibrium. The study of these quasiequilibrium configurations is 
justified in the view that the fully dynamical computations mentioned above are only in their infancy. In particular, 
they cannot follow more than a few orbits. Also they involve a rather coarse resolution of the stars, being performed 
in a single box with Cartesian coordinates. Another motivation for computing quasiequilibrium configurations is to 
provide valuable initial conditions for the dynamical evolutions 

The first quasiequilibrium configurations of binary neutron stars in general relativity have been obtained three 
years ago by Baumgarte et al. |2^,^], followed by Marronetti et al. pffl. However these computations considered 
synchronized binaries. This rotation state does not correspond to physical situations, since it has been shown that 
the gravitational-radiation driven evolution is too rapid for the viscous forces to synchronize the spin of each neutron 
star with the orbit [ p8|p9| ] as they do for ordinary stellar binaries. Rather, the viscosity is negligible and the fluid 
velocity circulation (with respect to some inertial frame) is conserved in these systems. Provided that the initial 
spins are not in the millisecond regime, this means that close binary configurations are well approximated by zero 
vorticity (i.e. irrotational) states. Irrotational configurations are more complicated to obtain because the fluid velocity 
does not vanish in the co-orbiting frame (as it does for synchronized binaries). We have successfully developed a 
numerical method to tackle this problem and presented the first quasiequilibrium configurations of irrotational binary 
neutron stars elsewhere pTJ . The numerical technique relies on a multi-domain spectral method pQ | within spherical 
coordinates. Since then, two other groups have obtained relativistic irrotational configurations: (i) Marronetti, 
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Mathews & Wilson I23J31I] b y m eans of single-domain finite difference method within Cartesian coordinates and (ii) 



Uryu, Eriguchi and Shibata [|23 24 by means of multi-domain finite difference method within spherical coordinates. 

The present article is devoted to the detailed presentation of our method, along with numerous tests of the numerical 
code, while the previous letter pj| gave only a sketch of the equations and some results about an evolutionary sequence 
built on a polytropic equation of state. In particular, that letter focuses on the evolution of the central density along 
the sequence in order to investigate the stability of each star against gravitational collapse. That study was motivated 
by the f 995 finding of Wilson et al. |(|[?J that the neutron stars may individually collapse into a black hole prior to 
merger. This unexpected result has been called into question by a number of authors (see Ref. |32| | for a summary of 
all the criticisms and some answers). Recently Flanagan ]33| has found an error in the analytical formulation used 
by Wilson et al. HQ. New numerical computations by Mathews and Wilson |34|) , using a corrected code, show a 
significantly reduced compression effect. 

The plan of the article is as follows. We start by presenting the equations governing binary stars in general relativity 
in Sect, [n] (hydrodynamics) and Sect. [II (gravitational field). The numerical method developed to integrate these 
equations is presented in Sect. IV. Section M is then devoted to the tests passed by the numerical code. Astrophysical 



results are then presented in Sect. VI for an evolutio nary sequence of irrotational binary stars constructed on polytropic 



equation of state of adiabatic index 7 = 2. Section VII contains the final discussion (comparison of our method with 
that used by other groups, conclusions about the tests) and future prospects. 

Throughout the present article, we use units of G = c = 1 where G and c denote the gravitational constant and 
speed of light. 

II. RELATIVISTIC EQUATIONS GOVERNING BINARIES IN CIRCULAR ORBITS 

Our treatment of binary neutron stars relies on the assumptions of (i) quasiequilibrium state (i.e. steady state in the 
co-orbiting frame), (ii) a specific velocity state for the fluid: either rigid or irrotational flow, (iii) the spatial 3-metric 
is almost conformally flat. In this section, we examine the assumptions (i) and (ii), without invoking assumption (iii), 



which will be introduced only in Sect. Ill 



A. Quasiequilibrium assumption 

In the late inspiral phase, before any orbital instability or merging of the two stars, the evolution of binary neutron 
stars can be approximated by a succession of circular orbits. Indeed when the separation between the centers of the 
two neutron stars is about 50 km (in harmonic coordinates) the time variation of the orbital period, P or b, computed 
at the 2nd Post-Newtonian (PN) order by means of the formulas established by Blanchet et al. |^5| is about 2%. The 
evolution at this stage can thus be still considered as a sequence of equilibrium configurations. Moreover the orbits 
are expected to be circular (vanishing eccentricity), as a consequence of the gravitational radiation reaction pq| . In 
terms of the spacetime geometry, we translate these assumptions by demanding that there exists a Killing vector field 
1 which is expressible as |l7j 

l = k+ftm, (1) 

where Q is a constant, to be identified with the orbital angular velocity with respect to a distant inertial observer, 
and k and m are two vector fields with the following properties: k is timelike at least far from the binary and 
is normalized so that far from the star it coincides with the 4-velocity of inertial observers with respect to which 
the total ADM 3-momentum of the system vanishes. On the other hand m is a spacelike vector field which has 
closed orbits and is zero on a two dimensional timelike surface, called the rotation axis, m is normalized so that 
V(m • m) ■ V(m • m)/(4m ■ m) tends to 1 on the rotation axis [this latter condition ensures that the parameter ip 
associated with m along its trajectories by m = d/dtp has the standard 2ir periodicity]. Let us call 1 the helicoidal 
Killing vector. We assume that 1 is a symmetry generator not only for the spacetime metric g but also for all the 
matter fields. In particular, 1 is tangent to the world tubes representing the surface of each star, hence its qualification 
of helicoidal (cf. Figure 1 of Ref. |17| ) . 

The approximation suggested above amounts to neglecting outgoing gravitational radiation. For non-axisymmetric 
systems — as binaries are — imposing 1 as an exact Killing vector leads to a spacetime which is not asymptotically 
flat fl37| . Thus, in solving for the gravitational field equations, a certain approximation has to be devised in order to 
avoid the divergence of some metric coefficients at infinity. For instance such an approximation could be the Wilson 
& Mathews scheme [ 
well as the trace of t 



that amounts to solving only for the Hamiltonian and momentum constraint equations, as 



re spatial part of the "dynamical" Einstein equations (cf. Sect. Ill A). This approximation has 
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been used in all the relativistic quasiequilibrium studies to date and is consistent with the existence of the helicoidal 
Killing vector field ([[]). Note also that since the gravitational radiation reaction shows up only at the 2.5-PN order, 
the helicoidal symmetry is exact up to the 2-PN order. 

Following the standard 3+1 formalism we introduce a foliation of spacetime by a family of spacelike hypersur- 
faces St such that at spatial infinity, the vector k introduced in Eq. ([!]) is normal to S( and the ADM 3-momentum 
in E t vanishes (i.e. the time t is the proper time of an asymptotic inertial observer at rest with respect to the binary 
system). Asymptotically, k = d/dt and m = d/dip, where ip is the azimuthal coordinate associated with the above 
asymptotic inertial observer, so that Eq. (Q) can be re-written as 

One can then introduce the shift vector B of co-orbiting coordinates by means of the orthogonal decomposition of 
1 with respect to the £< foliation: 

l = iVn-B, (3) 
where n is the unit future directed vector normal to S t , TV is called the lapse function and n • B = 0. 

B. Fluid motion 

We consider a perfect fluid, which constitutes an excellent approximation for neutron star matter. The matter 
stress-energy tensor is then 

T = (e+p)u®u+pg, (4) 

e being the fluid proper energy density, p the fluid pressure, u the fluid 4-velocity, and g the spacetime metric. A 
zero-temperature equation of state (EOS) is a very good approximation for neutron star matter. For such an EOS, 
the first law of thermodynamics gives rise to the following identity (Gibbs-Duhem relation): 

VP = \^ , (5) 



e + p h 
where h is the fluid specific enthalpy: 

, e + p 



rriBn 



(6) 



n being the fluid baryon number density and tob the mean baryon mass: tob = 1.66 x 10~ 27 kg. Note that for our 
zero-temperature EOS, h is equal to the fluid chemical potential. 

By means of the identity (jsj) , it is straightforward to show that the classical momentum-energy conservation equation 
V • T = is equivalent to the set of two equations |39],flO| : 

u • (V A w) = , (7) 

V • (nu) = . (8) 

In Eq. (0), w is the co-momentum 1-form 

w:=/iu (9) 

and V A w denotes the exterior derivative of w, i.e. the vorticity 2-form p9| . In terms of components, one has 

(V A w) Q(3 = V a wp - Vpw a = d a Wf3 - d/3W a . (10) 

The vorticity 2-form enters Cartan's identity which states that the Lie derivative of the 1-form w along the vector 
field 1 is 

Aw = l-(VAw) + V(l-w) . (11) 
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Because of the assumed helicoidal symmetry, £\\v = 0, so that Cartan's identity reduces to 

I- (V Aw) + V(l- w) =0 . (12) 

This equation reveals to be very useful in the following; this justifies the introduction of the vorticity 2-form. 
In particular, performing the scalar product of Eq. (|lj) by the fluid 4-velocity u leads to 

1- (V A w) • u + u- V(l- w) = . (13) 

The first term in the left-hand side vanishes by virtue of the equation of motion (]?]) , so that we obtain 

u-V(l-w)=0, (14) 

which means that the quantity 1 • w = h 1 • u is constant along each streamline. This is the relativistic generalization 
of the classical Bernoulli theorem. At this stage, it must be noticed that, in order for the constant to be uniform 
over the streamlines, i.e., to be a constant over spacetime, so that one gets a first integral for the fluid motion, 
some additional property of the flow must be required. In the following two sections, we explore two such additional 
properties: rigidity and irrotationality. 



C. Rigid rotation 

A rigid motion corresponds to synchronized stars (also called corotating stars). It is defined in relativity by the 
vanishing of the expansion tensor a p := (gj 1 + u a u ,1 )(g l3 l/ + U/3U I/ )V( l/ u /J ) of the 4-velocity u. In the presence of a 
Killing vector 1, this can be realized by requiring the colinearity of u and 1 : 

u = Al, (15) 

where A is a scalar field related to the norm of 1 by the normalization of the 4-velocity A = (—1 • l) -1 / 2 . Inserting 
relation (|l5|) into the equation of fluid motion (R) shows that the first term in Eq. ( |l2| ) vanishes identically, so that 
one gets the well known first integral of motion |4l| ] 

1 • w = const. (16) 

The second part of the equations of fluid motion, Eq. (^) (baryon number conservation), is trivially satisfied by the 
form ([D]) because 1 is a Killing vector. 



D. Irrotational flow 

As recalled in Sect. realistic binary neutron stars are not expected to be in synchronized rotation, but rather to 
have an irrotational motion. A relativistic irrotational flow is defined by the vanishing of the vorticity 2-form p9| : 

VAw = 0. (17) 

This is equivalent to the existence of a scalar field ^ such that 

w = W. (18) 

This is the relativistic definition of a potential flow . Note that the advantage of writing the equation for the fluid 
motion in the form rather than in the traditional form V • T = is that one can see immediately that a flow 

of the form (|l|) is a solution of (0). 



The second part of the equation of motion, Eq. (|8|), is satisfied by the potential flow ( 18) provided that ^ obeys to 
the equation 

-V- W + W- V( T ) =0 . (19) 
h \hJ 



Inserting the irrotationality condition (|17j) into Eq. (|12j) results in an equation showing the constancy of the scalar 
product 1 • w: 

1 w = const. (20) 

We therefore obtain the same first integral as in the rigid case (Eq. ( |l6| ) above). However note that the way to get it 
is different: no use of the equation of motion (0) has been made to obtain (20), contrary to the derivation of Eq. Jl€|). 



The first integral ( J16| ) f° r rigid motion has been known for a long time, at least since Boyer's work ]41| . To our 
knowledge, the version (p0h for an irrotational flow in presence of a Killing vector is due to Carter Gil 
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E. 3+1 decomposition 



The first integral (|l6|) , (pcj) , common to both the rigid and irrotational motion, is expressed in terms of the con- 
traction of a spacetime vector (1) with a spacetime 1-form (w). Going back to the 3+1 formalism mentioned in 
Sect. EI A, let us re-express it in terms of quantities relative to the hypersurfaces S t . Following Ref. we introduce 



the co-orbiting observer, whose 4-velocity v is the normalized symmetry generator: 

v = (TV 2 - B • B) -1 / 2 1 , (21) 

where the normalization factor has been deduced from Eq. (g). Note that in the rigid motion case, the co-orbiting 
observer and the fluid comoving observer coincide: u = v [cf. Eq. (p^)]. The 3+1 split of the 4-velocity v with respect 
to the Eulerian observer is 

v = r (n + U ), (22) 

where 

r = -n-v=(l-U -U )- 1 / 2 (23) 
is the Lorentz factor between the two observers and Uq is the orbital 3- velocity with respect to the Eulerian observer 



(n • Uq = 0). According to Eqs. (21) and (m, Uq is linked to the shift vector of co-orbiting coordinates by 



Uo = ~ • (24) 
Thanks to the second part of Eq. (|23|) , Eq. ( |2l| ) can be re- written as 

< 25 > 

The fluid motion can be described by the following orthogonal decompositions of u: 

u = T(v + V) =T n (n + U) , (26) 

where F = — v • u (resp. r n = n • u) is the Lorentz factor between the fluid and the co-orbiting (resp. Eulerian) 
observer, and V (resp. U) is the fluid 3-velocity with respect to the co-orbiting (resp. Eulerian) observer. In 
particular, v • V = 0, n • U = and 

U = -Lh-u, (27) 

-L n 



where 



h:=g + n®n (28) 



is the orthogonal projector onto the spatial hypersurfaces £*: h can also be viewed as the metric induced by g onto the 
hypersurfaces St. Performing the scalar product of Eq. (|22f ) with the second part of Eq. (^) leads to an expression 
of the Lorentz factor T in terms of quantities relative to the Eulerian observer only: 

r = r n r (i-u-u ) . (29) 

Similarly, performing the projection of the second part of Eq. (|2^) onto the hyperplane orthogonal to v results in the 
expression of the fluid 3-velocity V with respect to the co-orbiting observer in terms of the 3- velocities U and Uo, 
both defined with respect to the Eulerian observer: 

U [U -(U-Uo)n + U-Uo + (U-U )Uo-(Uo-U )U] . (30) 



1-U-U 



Note that in the case where U and Uo are aligned (U = Ue and Uo = Uoe, e being some unit vector in E t ) 
relation ( |30| ) reduces to the classical velocity- addition law of special relativity: V = (U — U )/(l — UUo)e', where 
e' = To(e + f/on) is the unit vector deduced from e by a boost of velocity Uq. In particular for U = Uo, which 
corresponds to synchronized binaries, V vanishes identically. 
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For irrotational binaries, U is related to the potential ^ by combining Eqs. (^), ( p"8| ) and ( pjj ) : 

U =i M > < 31 > 

where D is the covariant derivative associated with the metric h of spatial hypersurfaces St- Combined with the 
relation r n = (1 — U • U) -1 / 2 , this relation results in 

1 \ V 2 



— D*-D*) (32) 



We are now in position to write the 3+1 form of the first integral (16), (PU), common to both the rigid and irrotational 
motion. Substituting relation (^) for w and relation (|25| ) for 1 into Eq. (|16|) results in 

r 

hN — = const . (33) 
To 

We shall use actually the logarithm of this relation: 

H + v - In T + In T = const , (34) 



with the following definitions: 



and 



H:=\nh (35) 



u:=\nN . (36) 

These two quantities have immediate meaning at the Newtonian limit: H is the (non-relativistic) specific enthalpy 
and v is the Newtonian gravitational potential. The first integral of motion written as ( p4] ) coincides with Eq. (66) of 
Ref. [jf7|. The link with the alternative expressions derived by Teukolsky [[l9| and by Shibata [p0| for the irrotational 
case is performed in Appendix |A| N ote that lnT = for synchronized binaries, so that Eq. (|34|)simplifies somewhat. 
Note also that substituting Eq. (|29|) for T in Eq. (Q) leads to an alternative expression of the first integral of motion 
which contains only quantities relative to the Eulerian observer: 

H + v + lnT n + ln(l - U • U ) = const . (37) 

However in the following, we shall use only the form (|3^). 

Let us now turn to the 3+1 form the differential equation ([l9]) for the velocity potential \& of irrotational flows. 
Taking into account the helicoidal symmetry, Eq. ( |l9| ) becomes 

n D + Dn • = hT D U • Dn + n ■ D In A + U • Dr n \ + nhKT n , (38) 

where K is the trace of the extrinsic curvature tensor of the T, t hypersurfaces. This equation has been obtained 
by Teukolsky |l9| and independently by Shibata p0[ |. We refer to these authors for the details of the derivation of 
Eq. (||) from Eq. @. 



F. Newtonian limit 

At the Newtonian limit, the Eulerian observer is an inertial observer. Eqs. (^) and (||) show that B = — O^-, so 
that Eq. (^4|) for the velocity of the co-orbiting observer with respect to the inertial observer becomes 

U = !]xr, (39) 

where r denotes the position vector with respect to the center of mass of the system. The logarithm of the corre- 
sponding Lorentz factor tends to (minus) the centrifugal potential [cf. Eq. (p3|) 1 
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lnTo = x r) 2 . (40) 



The Newtonian limit of the first integral of motion (|34| ) for synchronized binaries (In T = 0) gives the classical 
expression 

H + v- -(CI x r) 2 = const , (41) 

where, as recalled above, H is the fluid specific enthalpy and v the Newtonian gravitational potential. 

In the irrotational case, the Newtonian limit results in the following fluid velocity with respect to the inertial frame 
[set h = 1 and T n = 1 in Eq. @] 

U = W , (42) 

where V denotes the standard 3-dimensional gradient operator, i.e. the Newtonian limit of the operator D introduced 
above. The first integral of motion (p^) reduces then to 



H + v+ ^(W) 2 - (£1 x r) • W = const. (43) 

We recognize the classical expression (compare e.g. with Eq. (12) of Ref. §]| or Eq. (11) of Ref. ||). The Newtonian 
limit of the continuity equation ( |38| ) reads 

nA* + Vn ■W = (J]xr)-Vn. (44) 
Here again, we recognize the classical expression (compare e.g. with Eq. (13) of Ref. |19| ). 

III. GRAVITATIONAL FIELD EQUATIONS 

A. A simplifying assumption: the conformally fiat 3-metric 

As a first step in the treatment of binary configurations in general relativity, a simplifying assumption can be 
introduced, in order to reduce the computational task, namely to take the 3-metric induced in the hypersurfaces E t 
to be conformally flat: 

h = A 2 f , (45) 

where A is some scalar field and f is a flat 3-metric. This assumption has been first introduced by Wilson & Mathews 
p8[ and has been employed in all the studies of quasiequilibrium relativistic binaries to date |2l]-[2(| . It has been also 
used in binary black hole initial data computations (see e.g. jf||-[l7|]). It is physically less justified than the assumption 
of quasiequilibrium discussed above. However, some possible justifications of ([HI ) are 

1. it is exact for spherically symmetric configurations; 



2. it is very accurate for axisymmetric rotating neutron stars 48 1; 

3. the 1-PN metric fits it; 

4. the 2.5-PN metric p9] deviates from it by only 2 % for two 1.4 Mq neutron stars as close as 30 km (in harmonic 



coordinates) |50 



B. Partial differential equations for the metric 

To benefit from the helicoidal symmetry, we use co-orbiting coordinates (t, x , x 2 , x 3 ), i.e. coordinates adapted to 
the Killing vector 1: d/dt = 1. Assuming the conformally flat form (|||) for h, the full spacetime metric takes then 
the forrrj^] 



1 Latin indices i, j, . . . run from 1 to 3. 
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ds 2 = -{ l N j - BiB l )dt 2 - 2Bidt dx l + A fydx 1 dx 3 . (46) 

We have thus five metric functions to determine: the lapse N, the conformal factor A and the three components B % of 
the shift vector B [see Eq. (H)]. Let us define auxiliary metric quantities: we have already introduced the logarithm 



of N, v, via Eq. (36); we introduce now the shift vector of non-rotating coordinates: 



and the quantity 



N = B + A , (47) 



:= ln(A/V). (48) 



At the Newtonian limit N = and (3 = 0. 

In the following, we choose the slicing of spacetime by the hypersurfaces S( to be maximal. This results in K = 0. 

The Killing equation V a i/3 + V pl a = 0, gives rise to a relation between the T, t extrinsic curvature tensor K and the 
shift vector N (via Eq. (0) and the relation Vn = — K — n <g> D In AT) 



K tj = — — (D*flj + D 3 B i ) = 1— \Vn 3 + V 3 N l - -f i3 V k N k \ , 

2N y ' 2A 2 N \ 3 J J ' 



(49) 



where V stands for the covariant derivative associated with the flat 3-metric f . Here and in the following, the index 

i of V is supposed to be raised with the metric f. Note that since d/dcj) is a Killing vector of the flat metric f, the 
second part of this equation stands also with N l replaced by B l . 

The trace of the spatial part of the Einstein equation, combined with the Hamiltonian constraint equation, result 
in the following two equations 

Au = 4nA 2 (E + S) + A 2 KijK i: > - VifV/3 , (50) 

A/3 = AttA 2 S + - l'-K,,K'' - \ (V^VV + ViPTp) , (51) 
4 2 \ / 

whereas the momentum constraint equation yields, by means of Eq. (|49|), 

AN 1 + \V (VjN j ) = -16irNA 2 (E +p)U l + 2NA 2 K 13 V 3 (3/3 - Av) . (52) 

In these equations, A := V V, is the Laplacian operator associated with the flat metric f , and E and S are respectively 
the matter energy density and trace of the stress tensor, both as measured by the Eulerian observer: 

E :=n-T-n = T 2 n (e+p)-p, (53) 

S :=h-T = 3p+(E + p)XJ -U . (54) 

The equations to be solved to get the metric coefficients are the elliptic equations (|50|)-(5^). Note that they represent 
only five of the ten Einstein equations. The remaining five Einstein equations are not used in this procedure. Moreover, 
some of these remaining equations are certainly violated, reflecting the fact that the conformally flat 3-metric (^) is 
an approximation to the exact metric generated by a binary system. 

At the Newtonian limit, Eqs. (jHl]) and ((5^) reduce to = 0. There remains only Eq. (0), which gives the usual 
Poisson equation for the gravitational potential v. 

C. Equations for the fluid with a conformally flat 3-metric 

In this section we explicitly write some equations for fluid quantities when the 3-metric takes the form (|45|). First 
the Lorentz factor ( p3| ) between the co-orbiting and Eulerian observers writes 

roHi-i 2 /,,* 1 ' 2 . (55) 



8 



For irrotational motion, the expression ( |32| ) for the Lorentz factor T n between the fluid and Eulerian observers 
becomes 

1 .._ _ \ 1/2 



ra= l 1 + ^ WV, ' $ ) ' (56) 



The corresponding fluid 3-velocity (J3l]) is 



The Lorentz factor T between the fluid and co-orbiting observer, which enters in the first integral of motion (|34|), is 
deduced from the above quantities via Eq. ( p9| ) : 

r = r n r„(i- A 2 n,u l u^ . (58) 

Let us now consider the continuity equation (jH). For a zero-temperature EOS, H can be considered as a function 
of the baryon density n solely, so that one can introduce the thcrmodynamical coefficient 

a inn 



The gradient of n which appears in Eq. ( p8[ ) can be then replaced by a gradient of H so that, using the metric (|4q), 
one obtains 

QHM> +V l HVit! = A 2 hT n U^iH + QH {V i ^V l {H ~ (3)+ A 2 hU^ l T 1 ^j . (60) 

The potential 'J is in fact dominated by a pure translational part. Therefore, we write, in each star, 

* =: *o + I, M ;,■>■' , (61) 
where Wq is the constant (translational) velocity field defined as the centra^] value of 

W l := A 2 hT n U l . (62) 

Then = V^o + and A* = A* , so that Eq. @ becomes 

r — % —i i — • — ( — W i — \ 

(HA}£ a + [(1 - C^)V H + (HV /3j Vi*o = {W % - W^ViH + (H (w ( )V l (H -(3) + — ViT n j . (63) 



The advantage to solve Eq. (p3|) instead of Eq. (60) is that the right-hand side of the former is much smaller than the 



right-hand side of the latter, due to the factor W — Wq, instead of W l , in front of ViH. 

D. Global quantities 

The total mass-energy content in a S t hypersurface is given by the Arnowitt-Deser-Misner (ADM) mass M, which 
is expressed by means of the surface integral at spatial infinity 

M = lib / fikfjl ( Wjhkl ~ Wkhjl ^ dSi (64) 



(see e.g. Eq. (20.9) of Rcf. |51|). In the case of the conformally flat 3-metric — A 2 fij, this integral can be written 



M = -77- i ^A 1 / 2 dS l . (65) 

27T ./„ 



The centers of the stars are denned in Sect. IV A. 
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By means of the Gauss-Ostragradsky formula, this expression can be converted into the volume integral of AA 1 / 2 . 
This last quantity can be expressed by subtracting Eq. (|5(]) from Eq. ( |5l| ) (recall that A = exp(/3 — v)"), so that 
Eq. ( |65| ) becomes an integral containing the matter energy density and the extrinsic curvature of St : 

M = J A 5/2 (e+ -l-KijK^ d 3 x . (66) 

Following Bowen & York |]52|| we define the total angular momentum in a St hypersurface as the surface integral 
at spatial infinity^ 

Ji = ^e ijk j {tPK* - x k K*) dSi = ^e ijk j> (x^A 5 K kl - x k A 5 K jl ) dS t , (67) 

where e^-fc is the 3-dimensional alternating tensor, x % are Cartesian coordinates, and the second equality follows from 
the fact that A = 1 at spatial infinity. As for M, this integral can be converted to a volume integral on S t . Using 
the momentum constraint equation D\K kl = 8n(E + p)U k and the fact that D\K kl = A~ 5 d(A 5 K kl )/dx l for the 
conformally flat 3-metric (]45|), one obtains the expression 

Ji = e ijk [ A 5 (E+p) x^U k d 3 x . (68) 



The baryon mass of each star is given by the integral on St of the baryon number density as measured by the 
Eulerian observer: — nu ■ n = T n n. In the case of the conformally flat 3-metric ([45[), this integral becomes 

A4 a> = m B / A 6 T n n d A x , a = 1, 2 . ( (>') ! 

star a 



IV. NUMERICAL METHOD 

The equations to be solved to get a relativistic binary system in quasiequilibrium are the elliptic equations (|50|)-(|52 
for the gravitational field, supplemented by the elliptic equation (J63J) for the velocity potential m the irrotational 
case. A cold matter equation of state, of the form 

n = n(H) e = e(H) p = p(H) , (70) 

must be supplied to close the system of equations. The thermodynamical quantity H has been privileged in the 



EOS setting ( ]T0| ) because it is that quantity which is involved in the first integral of motion (34). Altogether, these 
equations constitute a system of coupled non-linear partial differential equations. We solve this system by means of 
an iterative procedure. 



3 Note that contrary to the ADM mass, the total angular momentum hence defined is not asymptotically gauge invariant: it 
is defined merely as the 1/r 3 part of Kij within our coordinates; see York fc3| for a discussion. 
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FIG. 1. Coordinate systems used in the calculation. 



A. Coordinate systems and computational domains 

We use co-orbiting coordinates (t,X,Y,Z) of Cartesian type (i.e. — Sij), so that the line element ( fl6|) can be 
written 

ds 2 = -N 2 dt 2 + A 2 [(dX - B x dt) 2 + {dY - B Y dtf + (dZ - B z dtf] . (71) 

In these coordinates, the two stars have fixed locations and figures. Let us define the center of star no. a (a — 1, 2) as 
the location of the maximum enthalpy H (or equivalently maximum density e) in star a. Note that this center does 
not coincide with the center of mass of star a. We choose the coordinates (X, Y, Z) such that (i) the orbital plane is 



defined by Z = 0, (ii) the two stellar centers are located along the X axis and (iii) the rotation axis (cf. Sect. [I A) is 
located at X = 0, Y = 0. Let us then denote by Xn\ and X(q\ the X coordinates of the two stellar centers^. 

In order to describe properly the stellar interiors, we introduce two systems of Cartesian coordinates (xu\ , yi a \ , ) 
centered on the two stars by (see Fig. fy) 




X — Xq) \ x (2) - -{X- x {2)) 




= Y and { y {2) = —Y (72) 

= Z 

Note that the system (x^,y^,z^) is aligned with (X 7 Y,Z), whereas (2(2), U(2), z (2)) is anti-aligned (rotation of 
angle ir in the (X, Y) plane) with (X,Y,Z). This choice ensures that the companion of star no. a is located at 
X{ a ^ > for both stars. In particular, for equal mass stars, the descriptions of each star in terms of {x{ a ),y(a): z (a)) 
are identical. Furthermore we introduce spherical coordinates (r/ a \,6/ a \,ip/ a \) (a — 1,2) associated with each of the 
Cartesian coordinate systems (£( a )i J/(a)i z (a)) by means of the usual formulae. 

Since some of the equations to be solved are elliptic equations with non-compactly supported sources, the compu- 
tational domain must extend up to spatial infinity, i.e. must cover the full hypersurface S t , in order to put correct 
boundary conditions (flat spacetime). Any truncated computational domain ("box") would result in approximate 
boundary conditions, which inevitably would induce some error in the numerical solution. The technique to cover 
the full S t is to divide it in various domains, the outermost of it being compactified in order to deal with finite com- 
putational domain only pC| , Following the introduction of the two coordinate systems (?"( a ), 9(a)- f(a)) ( one centered 
on each star), we will actually use two sets of such domains: one centered on star no. 1, the other on star no. 2. 
The number of domains in each set is arbitrary, being simply equal or larger than 3. The list of the Nn\ + N/%\ 
computational domains is 



4 In all this article, indices or superscripts in () will label the two stars. 
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<n<i> -0(1) 77(1) 

U Q ■■■ U M {1) -1 ■■■ u N m -l 
T)<2> -r,(2> r) <2> 

where 

• Domain Pq (a = 1, 2) has the topology of a ball and contains the center of star a; it is designed thereafter as 
the nucleus. 

• M( a j (a = 1, 2) is the number of domains which cover the interior of star a. It obeys M( a \ > 1 and M( ) < iV< a \ —2. 
The outer boundary of domain 2?j^ x coincides exactly with the surface of star a. The topology of domains 
T>[ a \ . . . ,T>^y _ x is that of a spherical shell; these domains are designed thereafter as the shells. 

• Domains T>^g _ 2 cover the non-compactified part of the space outside star a; they are also called 
shells. The inner boundary of domain ^ coincides exactly with the surface of star a. 

• Domain Djy —1 is the most external one; it extends up to r = +00. We call this domain the compactified 

external domain (CED) since thanks to some compactification it will be mapped to a finite computational 
domain. 

Of course the two sets of domains overlap since U . . . U T>^ t _ x = T>^ U . . . U T>^ 2 _i = St. The various domains 



are represented in Fig. [2] for iVm = N, 



(2> 




FIG. 2. Domains used in the numerical computations, when Nn\ = N/2) = 3. The boundaries of domains T>^ , T>^ , T>^ 
and T) 1 ^ 1 are represented. Fhe outer boundaries of domains T>^ and T>^ are located at infinity and are therefore not plotted. 



Following the technique introduced previously [p0| , we define in each domain the computational coordinates (£, 9' , if') 
according tof] 



Lp = Lp 



(73) 



and 



in the nucleus: 



£ + (3£ 4 - 2e)F (9, Lp) + - (5f - 3£ 5 ) G (9, if) 



Z G [0, 1] 



(74) 



5 For the sake of clarity we omit here the star indices (a) on the spherical coordinates (r( a ), 6( a ), <P(a)) centered on star a. 
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in the shells (1 < I < N {a) - 2) 



r = ai 



i + \ (£ 3 - 3£ + 2) Fi{0, + \ (-e + 3£ + 2) G t (6, ip) 



■A, fe[-l,l]; (75) 



in the CED: 



2R C 



ED 



e £[-!,!]■ (76) 



In the above relations, ai and j3i are some constants, the functions Fi(6, ip) and Gi(6, ip) define the boundary of each 
domain: the outer boundary of the nucleus corresponds to £ = 1 and is given by the equation 

r = ao[l + F {0,<p)+Go(0,<p)] , (77) 

where F (9,ip) contains only odd Fourier harmonics in ip and Gq{9,<p) only even harmonics, the inner boundary of 
the shell no. I (1 < I < N/ a \ — 2) corresponds to £ = — 1 and is given by the equation 

r = a l [-l + F l (9,<p)]+p l , (78) 

whereas its outer boundary corresponds to £ = 1 and is given by the equation 

r = a l [l + G l (0,tp)]+0i. (79) 

Finally i?cED is the radius of the inner boundary of the CED, which is assumed to be spherical. 

B. Multi-domain spectral method 

In each domain, we expand the various physical fields in a series of basis functions of £, 9' and ip' . We use Chebyshev 
polynomials in £, trigonometrical polynomials or associated Legendre functions in 9' , and Fourier series in ip' . The 



by Nr (0 the number of coefficient in the £ expansion used in domain V\ a> , by Ng 1 ' (I) the number of coefficients in the 
0' expansion and by (I) the number of coefficients in the ip' expansion. We employ a collocation spectral method, 
which means that in each domain, a function can be described either by the coefficients of its spectral expansion or 
by its value at some particular grid points, called the collocation points p4j . The grids plotted in Fig. |^ show actually 
these collocation points. 

The spectral method amounts to reducing linear partial differential equations into a system of algebraic equations 
for the coefficients of the spectral expansions. We refer to |^5|,[56| for the details of this multi-domain spectral method 
and here simply recall some basic features: 

• As explained above, spherical-type coordinates (£, 9', ip') centered on each star are used: this ensures a much 
better description of the stars than by means of Cartesian coordinates. 

• These spherical-type coordinates are surface-fitted coordinates: i.e. the surface of each star lies at a constant 
value of the coordinate £ thanks to the mapping (£,9',(p') i— > (r,9,(p) defined by Eqs. (f74|)-(f75|). This ensures 
that the spectral method applied in each domain is free from any Gibbs phenomenon (spurious oscillations 
generated by discontinuities). 

• The outermost domain extends up to spatial infinity, thanks to the mapping (|76|). This enables us to put exact 
boundary conditions on the elliptic equations (|5(])-([32]) for the metric coefficients: spatial infinity is the only 
location where the metric is known in advance (Minkowski metric). 

• Thanks to the use of a spectral method in each domain, the numerical error is evanescent for analytical fields 
(e.g. density fields for a 7 = 2 equation of state), i.e. it decreases exponentially with the number of coefficients 
(or equivalently collocation grid points) used in the spectral expansions |p5|.p6[ . 



interested reader is referred to Sect. III. A of Ref. 1 30 1 for more details about these spectral expansions. Let us denote 
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C. Splitting of the metric quantities 



Having introduced two sets of computational domains (grids), we linearly split the metric potentials u, (3 and N l 
into 



V = V(i) + V( 2 ) = V{1) + V(2->1) = v (1^2) + ^<2> , 



(80) 



= P{l)+P{2) = + P{2->1) = P{1^2) + P{2) 



(81) 



(82) 



where the quantities labeled by "(a)" or "(6 — ► a)" (a = 1, 2, 6 = 3 — a) are defined at the collocation points of the 
domains T)j a ' centered on star a, and the quantities labeled by "(a)" and "(a — > 6)" represents the same physical 
field but described at different collocation points (those of domain sets T>\°^ and vf^ respectively), i.e. ^{1^2) = 
"(2^1) = V(2), etc... 

The basic idea underlying the splittings (p0[)-(p2|) is that for each metric potential, there are two primary quantities, 
those labeled by "(1)" and "(2)", which are "mostly generated" by respectively star 1 and star 2 and which we called 
the auto-potentials [the precise definitions are given by Eqs. (^3|)-(|85|) below]. The auto-potentials are obtained by 



solving the gravitational field equations, on domains for the "(1)" potentials, and on vf^ for the "(2)" ones. 
The quantities labeled by "(1 — > 2)" [resp. "(2 — > 1}"] are then merely representations of the "(1)" [resp. "(2)"] 
auto-potentials at the collocation points associated with the companion star. For this reason, we shall call them the 
comp-potentials. 

Following the splittings (|S0|)-(|82|), the gravitational field equations (|50|)-([52"|) are themselves split in two parts, 
and V(2) are thus defined as the solutions of the two equations 



Ai/ (o> = ^A 2 {E (a) + S {a) ) + Q (a) + Q {b ^ a) - WiU {a) [v'/3 (a) + (V l /3) (b ^ a> 
whereas [3^ and /3{ 2 ) are defined as the solutions of the two equations 



a = 1,2 (b = 3 - a) , (83) 



3 ^ .- 

A/?(a) = 4irA 2 S {a) + -{Q(a) + Q(b^a)) ~ ^i v (a) V ^(a) + (V ^) (b^a) 



a = 1,2 



3 -a) 



(84) 



and and are defined as the solutions of the two equations 



^4*) + 3 V ' ( V ^(a)) = -l^NA\E {a) + p {a) )U\ a) +Nkl ) (& [V,/3 (Q) + &j0){b-*a) 

-8 [V^ (a) + (V^) (b ^ a) ] ) , a = 1, 2 (b = 3 - a) . 



(85) 



In these equations, E^, S^, P{ a ): Ut a \ are the quantities relative to the fluid of star a only and defined respectively 



by Eqs. (£3J), (|5J), and (|27|). K 1 ' is defined from N\ a) according to 



(a) 



'(a) - 3. 



a = 1,2 , 



(86) 



so that the total extrinsic curvature is given by K %3 = {K l ,L + K Z ,L)/A 2 . Finally Q( a ) and Q(p^ a ) are defined by 



Q(a) ■= Aif^K^K^ , a 



= 1,2 



(87) 



Q(b-~) ■= A'fikfjtK^KH^ , a = 1,2 (b = 3 - a) 



where K?^, is the same physical field than K % {L but numerically described at the collocation points of the domains 
X> ; , K % A being given at the collocation points of the domains Dj . 
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It is straightforward to check that adding the two equations (83) results in Eq. (J50|), adding the two equations fej ) 
results in Eq. ( |5l|) and adding the two equations ( [S5"l) results in Eq. (^). Therefore, having obtained solutions V{ a ), 
f3( a ) and of the equations (|83|)-(|85|), we can form the solution of the gravitational field equations (|50|), ( [sil ) and 
© via Eqs. 

The advantage of solving the system of 2 x 5 = 10 PDEs (||)-@, instead of solving the system of 5 PDEs (|5(j)-([52]), 
is that the source terms (right-hand-side) of the former are most ly con centrated on one of the two stars and therefore 
well described by one of the two domain sets introduced in Sect. IV A. This is not true for the source terms involving 
the comp-potentials "(& — > a)". However these terms enter only via quadratic combinations in which each of them 
is multiplied by the gradient of an auto-potential term, which is small where the comp-potential is large, so that the 
product of the two is smaller than the other sources terms, such as the scalar product of gradients of auto-potentials. 
The same considerations hold for Qi^^a) which appears to be much smaller than Q( a y According to these remarks, 



Eq. ( |S3| ) for is naturally solved on domains "D\ ; , Eq. (|S3|) for 1/(2) is solved on domains T>^' , and more generally, 
each equation for an auto-potential is solved onto the domains set centered on the corresponding star. 

Once the auto-potentials are known (at a given step of the iterative procedure described in the next section), there 
remains to compute the corresponding comp-potentials. This means that given e.g. vi^ at the collocation points of 

domains T>^\ one has to compute its values V(i^2) at the collocation points of domains T)^. One may think first 
to use some interpolation technique since the two sets of domains overlap. But this will necessarily introduce some 
"numerical noise" . We will proceed differently, taking advantage of the use of a spectral method. Indeed, the values 
of the field at the collocation points of domains T)\^ is not the only numerical representation of un\ we have 
at our disposal. We can use the alternative representation by the set of coefficients of its spectral expansion in each 



(2) 



domain (0 < I < Nn\ — 1) (cf. Sect. IV B ). By means of this spectral expansion, we can compute the value of 



z^i) at any point in the domain T>\ X \ not necessarily a collocation point. Hence, given a collocation point (£i,9j, <p' k ) 

of domain , we first compute the corresponding physical spherical coordinates (v"(2)) ^(2)1 ¥(2)) vm Eqs. (|74])-(|76|), 
then the corresponding Cartesian coordinates (^(2) , 2/{2)j z (2))] these latter are translated into Cartesian coordinates 
(xn\,yh\, Zn\) via Eq. (|72|). We finally obtain the corresponding spherical coordinates (»~{i),0(i),</?(i)) centered on 

star 1. We then determine in which domain T>^ this point is localized and to which value of the coordinate £ it 
corresponds by inverting the relations (|74|)-(|76|). Then we may use the spectral expansion of to get the searched 
value: 



^(l-2)(fo,£i)0i)¥>fc) 



E 

k=0 



-w fl (1> (Q- 

E 

.7=0 



E vikjiXkjiiCj I 9jy(0(i>) 



*fc(v<i>) . 



(89) 



where the vikji are the coefficients of i>/i\ in domain X*; 1 ^, Xkji denotes the basis functions in £ (typically Chebyshev 
polynomials), Qkj the basis functions in 9 (typically cosn9 or smn0) and $fc the basis functions in <p (Fourier series). 
These functions depend on the type of domain (nucleus, shell or CED) and are described in details in Sect. III. A of 
Ref. M. 



D. Iterative procedure 

Within our procedure, a quasiequilibrium binary neutron star configuration is obtained by specifying: 
1. the equation of state rt7G) for each star; 



2. the rotation state: either rigidly rotating (synchronized binaries, Sect. 1IC) or irrotational flow (Sect. II D ) 



3. the coordinate distance d :— \Xi2\ — -X"(i)| between the two stellar centers; 

4. the central enthalpies H?^ and in each star, or equivalently, via (f7(j|), the central density in each star (with 
our definition of the stellar center, this coincides with the maximum density) . 



As we discuss below, item can be replaced by the specification of the baryon mass of each star. 
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1. Initial conditions 



The above parameters being set, we start by computing initial conditions for the iterative procedure. These 
initial conditions are constituted by two numerical solutions for spherically symmetric static isolated neutron stars, of 
respective central enthalpy and H? 2 y Mm an d M^) being the gravitational masses of these spherical symmetric 
models, we set the X coordinates of the two stellar centers according to the Newtonian-like formulas: 



X 



M, 



(!) 



2) 



These coordinates will remain fixed during the iteration. Only the location of the rotation axis X rot , initially set to 0, 
will change (see Fig. [j]). Accordingly the formulae ( |90| ) have no physical meaning whatsoever. They can be viewed as 
the setting of the origin of the coordinate system (X, Y, Z) . The location of this origin is a priori arbitrary, only the 
distance d between the two stellar centers having a physical meaning; the setting (BQ) simply insures that this origin 
is not too far from the rotation axis. 

The angular velocity Q is initialized according to a formula for second order Post-Newtonian spherical stars fl35|,p7f : 



and 



X 



(2) 



M (1> 



Mm + M, 



(90) 



O 2 . = 

d 3 



1 - 



+ 



d 

Mini 



2R 2 



■7' 



d 2 

69 11 R 2 



12 R 4 2 
25 ~¥ 7 



17 R 4 



4 d 2 7 + 25 d 4 7 



(91) 



where M- m \ := Mm + M( 2 ), R is the coordinate radius of one of the two stars^J (which is spherical initially) and 
7 = 7irrot := for irrotational binaries, whereas 7 = 7 CO rot := 5I( a ) /{ZM^R 2 ^) for corotating binaries, I( a ) being 
the moment of inertia of star a. For this last quantity, we use as an ansatz the exact value for a Newtonian n = 1 
polytrope, which results in 7 coro t = 5/3 (1 — 6/-7T 2 ), independent of a. 

The metric auto-potentials are initialized as follows: and /3( a ) are set to the values of v and [3 for the static 
spherical models. The shift iV"? . is initialized to the first-order Post-Newtonian value for spherical incompressible 
binaries (this value can be obtained by taking the limit for a spherical star of the equations presented in Ref. [0]) 



N 



7 



(a) 



w, 



(a) 



V X(a> + V Wf a 



(a)" 



1,2 



with 



6A/ (a> O ini !i 
£ < a > (l+Af (a> /M (b> )fl <a) 

4M (a> n in id 
£ < a > (l+M (a> /Af <b> )r (a> 



3K- 



(e(i) := -1, £( 2 ) := 1, b = 3 - a) and 



X(a) 



(l+M (a) /M( b) )H (a> 
4M {a) n ini dR 2 (a) y {a) 
5(1 + M {a) /M {b} ) ^ 



V{a) 



3r U 



for 
for 

for 
for 



r {a ) 
r {a ) 



> R(a) 



(92) 



(93) 



-<a> 



<a> 



(94) 



We refer to Sect. IV A for the definition of the coordinates X, Y, Z, and involved in these formulas. It 



appeared that the initial shift vector given above results in a too large angular velocity in the first steps. Therefore, 
we artificially lower it by multiplying it by 0.6. From this initial value of the shift, we get initial values of KrL vm 
Eq. (pq), and initial values of B and Uo via Eqs. ( ^7|) and (M). 

Regarding the fluid quantities, the 3-velocity U is initialized to Uo in the s ync hronized case, whereas in the 
irrotational case, ^0 is initialized to zero and \& is initialized accordingly via Eq. (J6l|); the Lorentz factor T n is then 
initialized via Eq. ([36]) and the 3-velocity U is initialized according to Eq. <E%. We get then initial values of the 
Eulerian energy density E and the trace of stress tensor S via Eqs. (|53| ) and (gjj). In these equations, we use for the 
proper energy density e and pressure p the values of the spherically symmetric initial stellar models. 



Equation ( |9l| ) is valid only for equal-mass stars binaries. There also exists a more complicated formula for stars with different 
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2. Description of one step 



At a given step, we start by determining the value of the orbital angular velocity 0, and the value of the X coordinate 
of the rotation axis, X mt (see Fig. |l|), by taking the gradient along X of the first integral of motion (|34|). Demanding 
that the enthalpy H be maximal at the center of each star (our definition of center) , this results in the two equations 



d 



dX 



hiTo 



d 



(*<o>,0,0) 



dX 



(^ + inr) 



(*<a>.0,0) 



o=l,2 , 

where lnTp can be expressed in terms of f2 and X rot thanks to Eqs. ([2 3D, ([24}), ( |45| ) and (f 

"(or + a^) 2 + (npr - x^) - n y ) 2 + {n z ) 2 



lnTn 



--Mi 



A 2 



(95) 



(96) 



Inserting this relation into Eq. ( [95] ) and setting Y = Z = 0, X = X^ or X( 2 ) results in a system of two equations 
for the two unknowns f2 and X Iot . This system is solved by standard methods. Having determined f2 and A rot , we 
can compute the components of the orbiting velocity Uq, via Eqs. (p3) and (Et|) : 



1 

A 



(n(x - x mt ) - n y ) 



' N 



(97) 



where A, N x , N Y and N z are the values of the lapse function and the components of the non-rotating-coordinates 
shift vector taken from the previous step. From Uo, we of course compute the Lorentz factor Tq by Eq. ([23]). The fluid 
3-velocity with respect to the Eulerian observer, U is set to Uo in the synchronized case, where in the irrotational 
case, \& is deduced from Uo and the previous step value of 4"o via Eq- (|6l|); the Lorentz factor T n is then computed 
via Eq. (E3J3J) and the 3-velocity U follows from Eq. (pTl). The Lorentz factor T between the fluid and the co-orbiting 



observer is deduced from the above quantities via Eq. (58). 

The elliptic equation ([33]) for ^0 is then solved by the numerical method described in Appendix |b} 
We then adapt the computational domains to the stars as follows. The first integral of motion (E 

following the splitting rtgQ) 



H = H (a) + V ta) + $ <a),oxt - v (a) ~ $(a>,ext , 



is written, 



(98) 



where we have introduced the "external" potential 



(a) ,ext 



= ^(b^a) - hiL + lnT 



(99) 



and the superscript "c" stands for values at the center of the star. First, we rescale the auto-potential vi a \ by a factor 
a 2 to make sure that the enthalpy vanishes at the point 0^ = w/2, (p^ — on the external boundary of domain 



V 



(a) 



M 



H' 



(a) 



(a) ,ext 



(a) ,ext 



(a) (a) 



(100) 



where the superscript "s" stands for values at the point £ = 1, 6/ a \ — 7r/2, ipi a \ — of domain T> 



t } ,. When the 

iteration converges, a tends to 1. We then replace by a 2 v^ in Eq. (|9^) to get the enthalpy field in all space. 
Following the technique described in Ref. |3Q ], we then compute new functions Fi(6, if) and Gi(9, ip) in the mappings 
(Q) and ([75|) in order to make the outer boundary of domain T>^ x coincide exactly with the surface of the star. 
Since the collocation points of the new mapping do not coincide (in the physical space, described by the coordinates 
( r (a)i ${a) i f{a))) with that of the previous mapping, the values of the enthalpy field at the new collocation points have 
to be computed. The details of this computations are explained in Sect. V.A of Ref. |50| 

From this new value of H , wc compute the fluid proper baryon density n, proper energy density e and pressure p 
via the EOS (|70|). We then get the Eulerian energy density E and the trace of stress tensor S via Eqs. (53) and (|54|). 
These last quantities are subsequently used to evaluate the source terms of the elliptic equations (|8^)-(85) for the 
gravitational potentials. These equations are solved by means of the multi-domain scalar and vector Poisson solvers 
for non-compact sources described in details in Refs. [BftB6|. In particular the vector Poisson equation (p9) for the 
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auto shift N 1 ,-. is reduced to a set of 4 scalar Poisson equations according to the scheme used by Shibata and Oohara 




Before the beginning of a new step, some relaxation is performed onto the enthalpy field and the auto-potentials, 
according to 

Q J <- XQ J + (1 - X)Q J - 1 , (101) 

where Q stands for any of the fields H, V( a )i P(a) ancl ^\ a ) ( a — 1)2), J (resp. J — 1) labels the current step (resp. 
previous step), and A is the relaxation factor, typically chosen to be 0.5 for H and 0.65 for the auto-potentials. 

For appreciably relativistic configurations, it appeared that the above relaxation is not sufficient to ensure the 
convergence. In this case, we update the comp-potentials not every step but every m steps, with typically m = 8. 
This slows the convergence but enforces it. 



3. Convergence to a given baryon mass 



In order to compute evolutionary sequences of binary neutron stars, one should be able to compute configurations 
for a given baryon mass, since this quantity is conserved during the gravitational-radiation driven evolution of the 
system. The baryon mass, given by Eq. ((6^), is not a natural parameter we can set in our procedure. As stated above, 
the freely specifiable parameters which fix one configuration are the coordinate distance d between the two stellar 
centers and the central enthalpies and -H^ 2 ) m eacn star. However, we can use the iteration procedure itself to 
make the final configuration have a specified baryon mass. Indeed, since the baryon mass is an increasing function 
of the central enthalpy (at least for the stable stars we are studying), we multiply at each step, the central enthalpy 
H%\ by the factor 



2< 

where £ is the relative discrepancy between the actual baryon mass at the considered step, J and the wanted 

baryon mass : £ := Mg ^/Mg — 1. When the iterative procedure converge, the factor rj tends to one. The 
same treatment is performed on star 2. 



4- End of the iteration 



To control the convergence of the iterative procedure, we introduce the relative difference between the enthalpy 
fields of two successive steps: 

\H J (xi) -H J - 1 (x i )\ , s 

where the summation is extended to all the collocation points X{ inside the star and J is the step label. 

We use typically SH = 10 -7 as a criterion to end the iteration. For very high precision calculations (check with 
analytical solutions, see below), we use instead SH = 10 -12 . 



E. Treatment of cusps 



For very close configurations, an angular point (cusp) may appear at the surface of the stars, similar to that in 
the Roche lobe at the Lagrange point L\ in the Roche problem. At this point the enthalpy gradient dH/dr vanishes 
in the direction of the companion. The surface of the star is then no longer smooth and surface cannot be decribed 
by the differentiable mapping (f74|)-(|75|), because the functions Fi(6, ip) and Gi(9, ip) are assumed to be expandable in 
cos(fc#) or sin(fc6*) series and in Fourier series in tp, which implies that they are smooth functions of (9, tp). 

The solution to this problem consists in freezing the adaptation of the mapping to the stellar surface when the 
enthalpy gradient becomes too small at the surface point which faces the companion star. More precisely, we define 
the ratio 
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x : = ( S°r mp ' (io4) 

where the index "eq, comp" stands for the value at the point (0{ a ) — tt/2, (p/ a \ = 0) on the stellar surface, whereas 
the "index" pole stands for the value at the point (8^ — 0, ip/ a \ = 0) on the stellar surface. When x passes below a 
certain threshold Xfr during the iteration process, we stop the adaptation of the mapping to the surface of the star. 
Xfr is chosen typically chosen between 0.3 and 0.55. 

In this case, a Gibbs phenomenon is present. The accuracy of the calculation is then lower than when the mapping 
is adapted to the surface of the star. However, since the difference between the stellar surface at the domain boundary 
is pretty small, the Gibbs phenomenon is rather limited. 

For irrotational configurations, the non-coincidence of the stellar surface with a domain boundary introduces a 



small error in the resolution of the equation (63) for the velocity potential \l/o by means of the technique explained in 
Appendix |B]. 



F. Numerical implementation 

The numerical code implementing the method described above is written in the LORENE (Langage Objet POUR 
LA RElativite NumeriquE) language | p9[ , which is a C++ based language for numerical relativity. A typical run 
makes use of 6 domains (N/-n — N{ 2 ) — 3 and Mm\ — M/2) — 1), with N r x Ng x N v — 33 x 21 x 20 coefficients in 
each domain. The corresponding memory requirement is 232 MB for an irrotational configuration. A computation 
involves ~ 250 steps, which takes 14 h on one CPU of a SGI Origin200 computer (MIPS R10000 processor at 180 
MHz). If the number of coefficients is lowered to N r x Ng x N v = 25 x 17 x 16, the memory requirement and CPU 
times becomes respectively 100 MB and 6 h 30 min. 

Note that due to the rather small memory requirement, runs can be performed in parallel on a multi-processor 
platform. This is especially useful to compute sequences of configurations. 

Both Newtonian and relativistic configurations, either corotating or irrotational, are calculated by the same code. 
Only the parts of the computation specific to one of these four cases are treated by different branches of the code. 



V. TESTS OF THE NUMERICAL CODE 

After constructing a numerical code for calculation of binary neutron stars, we must assert its validity by performing 
self-consistency checks and comparing the results with those of analytic solutions or those of previous numerical works. 
The plan of the tests of the numerical code is as follows. For the irrotational configurations: 

1. Check the convergence of the iterative procedure. 

2. Check the convergence of the global quantities when increasing the number of coefficients of the spectral method. 

3. Check the decay of the relative error in the virial theorem for Newtonian binary systems when increasing the 
number of coefficients of the spectral method. 

4. Check the agreement with some analytic solutions for Newtonian binary systems. 

5. Check the agreement with previous numerical solutions for Newtonian binary systems. 

6. Check the coincidence of the results of the purely Newtonian calculation with those of general relativistic one 
with small compactness. 

and for the corotating configurations: 

1. Check the agreement with previous numerical solutions of relativistic binary systems. 

For the purpose of the test computations, we consider identical star binary systems with the polytropic equation of 
state 

1/(7-1) 



n(H) 



(exp(iJ) - 1) 



(105) 

7 K 

p(H) = nniH) 1 (106) 

e(H) = — ^-rn(Hy + m B n(H) , (107) 
7-1 

where 7, k and me are some constants. For me we will use ma = 1.66 x 10 -27 kg (mean baryon mass). 
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A. Self-consistency checks 



First of all, we show in Fig. || the convergence of the iterative procedure described in Sect. IV E . This convergence 
is meas ured by means of the relative difference 8H between two successive steps values of the enthalpy field, as given 
by Eq. ( 103| ). The bump around the 70th step co rresponds to the switch on the procedure of convergence towards 



a given baryon mass, as described in Sect. IV D 3. One can notice systematic oscillations in the convergence curve 
every 8 step s. They result from the fact that the comp-potentials are updated only every 8 steps, as discussed in 
Sect. [VD2. We stop the iterations when the convergence has reached the 8H = 10~ 7 level (dashed horizontal line 
in Fig~|^ 
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FIG. 3. Convergence (measured by the relative difference SH in the enthalpy field between two successive steps) of the 
iterative procedure for one of the irrotational models with N r x No x N v — 33 x 21 x 20 collocation points. The bump around 
the 70th step corresponds to the switch on of the procedure of convergence towards a given baryon mass. 



Next, we show the convergences of the global quantities (i.e. ADM mass and total angular momentum) f or one 
configuration when the number of spectral coefficients (or equivalently the number of collocation points, cf. Sect. IV B) 
is increased. Furthermore, we present the convergence of the relative change in central energy density along a 
sequence when we increase the number of spectral coefficients. The calculations are performed for the case of 7 = 2, 
k = 0.0332 Pnuc - 2 (Pnuc : = 1-66 x 10 17 kg m~ 3 ); the baryon mass is Mb = 1.625M Q , which corresponds to the 
compactness M/R = 0.14 for isolated spherical stars. The coordinate separ ation d is taken to be 60 km. Six domains 



are used, with the following parameters (using the notations of Sect. |VAj): Nm = Np) — 3, M(i) = M, 



the same number of coefficients in each domains: Nr (0) = iVr X '(l) 



N r ,N^(0)=N^(l) 



'<2> 



: 1, with 
No and 



n£ } {0) 



N v . 



The ADM mass and total angular momentum are shown in Figs. [| and || as functions of N r . We used the following 
numbers of spectral coefficients: N r x Ng x N v = 9 x 7 x 6, 13 x 9 x 8, 17 x 13 x 12, 25 x 17 x 16, 33 x 21 x 20 and 
41 x 25 x 24. In Fig. §, we give the relative change in central energy density along a quasiequilibrium sequence for 
various numbers of spectral coefficients: (N r ,Ng,N v ) = (9,7,6), (13,9,8), (17,13,12), (25,17,16), (33,21,20), and 
(33, 25, 24). We find that these global quantities settle to a constant value (variations below ~ 10~ 5 ) for N r > 25. 
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Convergence of ADM mass 
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FIG. 4. Convergence of the ADM mass for one of the irrotational relativistic models, as the number of collocation points 
(or equivalently of spectral coefficients) is increased. 
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Convergence of total angular momentum 
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FIG. 5. Same as Fig. H but for the total angular momentum. 
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Convergence of relative change in e c 
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FIG. 6. Convergence of the evolution of the central energy density along a quasiequilibrium sequence, as the number 
7V r x No x N v of collocation points (or equivalently of spectral coefficients) is increased. 



B. Tests in the Newtonian regime 



In order to test Newtonian calculations, we compute a Mb = 10 3 M Q Newtonian sequence based on a polytropic 



equation of state with 7 = 2 and n = 0.0332 p n 



In this case, the central baryon density and the radius of infinitely 



separated stars be com es 1.08 1 x 10 p nuc and 20.57 km, respectively. Note that the Newtonian limit of the polytropic 



equation of state ( |105| )- ( 107 ) is obtained for H <C 1 and reads 



n(H) = 

p(H) = 
e(H) = 



7 - 1 m B 

7 K 

tub n{H) . 



H 



1/(7-1) 



(108) 

(109) 
(110) 



1. Virial theorem 



A useful method to estimate the global numerical error in Newtonian computations is to calculate the relative error 
in the virial theorem. This latter states that 2T + W + 3P = 0, where T, W, and P denote respectively the kinetic 
energy of the binary system, its gravitational potential energy and the volume integral of the fluid pressure. We 
therefore define the virial error as 

\2T + W + 3P\ , . 

Error = ^ • (HI) 
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Convergence of relative error in virial theorem 

Newtonian hrotational calculation: y=2, K=0.0332 p V, M B =0.001 M , 
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FIG. 7. Relative error in the virial theorem along an evolutionary sequence, for various numbers N r x Ne x N v of collocation 
points (or equivalently of spectral coefficients). 

This error estimator is shown along a constant baryon number sequence in Fig. |^. In order to check the convergence 
of the numerical method, we present various cases of increasing number of spectral coefficients: N r xNgX N v = 9x7x6, 
13 x 9 x 8, 17 x 13 x 12, 25 x 17 x 16, 33 x 21 x 20 and 33 x 25 x 24. We used SH = 1CT 12 as the criterion to end 
the iteration. It is found from Fig. [7] that for large separations the relative error converges to 10~ 12 when the number 
of spectral coefficients is increased, which is of the order of SH. Anyway, one cannot go much further, even if SH 
is lowered significantly, because of the use of 15 digits numbers (double precision) and the resulting accumulation of 
round-off errors in the arithmetical operations. Besides, we notice in Fig. |?] the appearance of a rapidly increasing 
error for closer separations. Note however that at the point of closest approach (cusp point) this error is below 10~ 4 
(except for the low numbers of spectral coefficients), which is very satisfactory. Finally, we see on Fig. ^ that when 
we increase the number of polar and azimuthal collocation points fixing the number of radial ones, the relative error 
in the virial theorem becomes better around intermediate separations. 



2. Comparison with analytic solutions 



Until recently the only analytic solutions for binary stars were constructed with incompressible fluid and belong to 
the so-called families of Roche- Riemann or Darwin- Riemann ellipsoids 7 |50 61 1 (see Ref . |6^] for a good introduction 
to these ellipsoidal solutions). We have presented elsewhere ]3(J the comparison with Roche ellipsoids (the sub-class of 
Roche- Riemann ellipsoids constituted by synchronized systems) , as a validation of our multi-domain spectral approach 
with surface-fitted coordinates. As can be seen from Fig. 6 of Ref. | p0| the numerical error is decreasing exponentially 
with the number of spectral coefficients (the so-called evanescent error typical of spectral methods), reaching 1CP 9 
for a Roche ellipsoid with axis ratios a^/ai = 0.75 and 03/01 = 0.68. 

The case of compressible fluid bodies has been investigated recently by Taniguchi & Nakamura , who have 

obtained semi-analytic solutions for equilibrium sequences of irrotational binary polytropic stars in Newtonian gravity. 
For an equal-mass star binary system with 7 = 2, they have produced the following simple equations for the total 
energy E, the total angular momentum J, the orbital angular velocity and the relative change in central baryon 
density: 



E = 



GM 2 
Ro 





'Ro\ 






+2 ( 




v d J 





15 



' Rn 



j = -Md z n 



1 + higher term than 0[ y — 



(112) 
(113) 



7 Note however that these solutions are not exact for the gravitational potential of the companion must be truncated to the 
second order to get perfectelly ellipsoidal shapes 
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o 2 = 



5p c 



2GM 

Pc - PcO 
PcO 



1 + 6 



3 "IKS)' 
2tt 2 W / ' 



(114) 
(115) 



where d is the separation, p c the central density, R$ the radius of the spherical star of mass M (i.e. the radius at 
infinite separation) and p cQ the central density of this spherical star. These equations are exact up to 0((Ro/d) e ) an d 
are very valuable to check the validity of the Newtonian limit of our code, in particular the solution of equation ( |6^ ) 
for the velocity potential. Indeed the Darwin-Riemann solutions could not have been used for testing this important 
part of the code because Eq. ( |63| ) is degenerate for an incompressible fluid (£ = oo). 

First, we compare our numerical results with Taniguchi & Nakamura's analytic solutions for global quantities (such 
that total energy, total angular momentum, orbital angular velocity, and relative change in the central baryon density) 
along an evolutionary sequence in Figs. || - [ll]. For the numerical computation, we use N r x Ng x N v = 33 x 25 x 24 
spectral coefficients in each domain and the criterion 5H = 10~ 12 to end the iterations. It is found from these figures 
that the numerical results agree very well with the analytic ones. Note that the analytic solution ends at the contact 
point, whereas the numerical one ends before (when a cusp appears at the stellar surface, cf. Sect. IV E). However the 
analytic solution, based on an expansion up to 0((Ro/d) 6 ), loses its accuracy for very close separations and cannot 
be used to test the code in this regime. 

Total energy compared with analytic solution 
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FIG. 8. Total energy compared with Taniguchi & Nakamura's analytic solution j6^,Q along an evolutionary sequence. 
Solid and dashed lines denote the results of numerical and analytic calculations, respectively. 



Total angular momentum compared with analytic solution 
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FIG. 9. Same as Fig. M but the total angular momentum. 



Orbital angular velocity compared with analytic solution 
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FIG. 10. Same as Fig. q but the orbital angular velocity. 



Relative change in central baryon density 
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FIG. 11. Same as Fig. |s| but the relative change in central baryon density. 

In order to investigate the discrepancy between the results from the numerical code and those from Taniguchi & 
Nakamura's analytic solution, we present the relative differences on global quantities as functions of the separation in 
a log-log plot in Fig. |l2|. The relative differences are defined as follows: 



GM 2 /R ' 

'J num "ana 

Md 2 n K e P /2' 

^num ^ana 
^Kcp 

\Sp 

c:num 

-S P 

c:ana 1 5 



(116) 

(117) 

(118) 
(119) 



where f^Kcp is the Keplerian velocity for point mass particles: 



25 



Kcp 



/2GAf\i/2 



(120) 



Two reference lines, proportional to (d/Ro) and (d/Ro) , have been drawn in Fig. [12| in order to check the slope 
of the results easily. 

It is found that for separations closer than d/Ro = 10, the discrepancies between numerical and analytic solutions 
for the energy E and the relative change in central density 5p c are both proportional to ~ (d/Ro)~ 9 , and become 
~ (d/i?o)~ 13 around d/Ro ~ 3. At first glance, this agreement betwe en th e nu meric al and analytical solutions seems 
too good, because we know that the next order term missing in Eqs. flll2| ) and (115) is 0(~ (d/ R )~ 8 ). We interpret 
the fact that this term does not show up in Fig. [l2] by the fact that it is produced by the octupole deformation, which 
should be very small. Of course, for separations much larger than d/Ro = 10, the term proportional to ~ (d/Ro)~ 8 
would dominate the inclination of the lines. 

For the residual terms of the angular momentum J and the orbital angular velocity f2, we can see that while they 
are proportional to ~ (d/Ro)~ 7 around d/Ro ~ 10, the term proportional to ~ (d/Ro)~ 8 dominates for separations 
closer than d/Ro ~ 6, and finally goes up to ~ (d/R )~ 12 . We can explain this dependence as follows. First, the high 
order expansion of Q can be written as 



n = n 



Kcp 



1 



■-no 



(121) 



Note here that the second ter m in side the brackets comes from the quadrup olc d eformation of the star and is included 
in the analytic solution [Eq. (114)]. After subtracting the analytic solution (114) from expression (121), there remains 
the term 0((d/ Ro)~ 7 ) as a leading one. Therefore it dominates the behavior of the curve of the relative difference in 
around d/Ro ~ 10. 

On the other hand, the angular momentum is expandable as 



2 



1 



(122) 



This means that after subtracting the analytic solution (113) and normalizing by Md 2 £lKep/2, the leading term of J 
comes from fi, because this latter has a term proportional to (d/Ro) . This explains why the discrepancy curves for 
J and O have almost the same behavior. 

From the above discussion about the slopes of Fig. [jj] curves, we can conclude that the numerical solution agrees 
with the semi-analytical one within the accuracy of this latter, i.e. the in crea se of t he discrepancy when the separation 
decreases is due to missing (high order) terms in the analytic solution (112)-(115). Finally, we see from Fig. [l^ that 
even if the baryon mass is changed by a factor larger than 10 3 , the numerical and analytical solutions agree in the 
very same manner. 
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FIG. 12. Relative differences in total energy E, total angular momentum J, orbital angular velocity Q, and relative change 
in central baryon density Sp c when comparing the numerical solution with Taniguchi & Nakamura's analytic solution [p3lp4] 
along an equilibrium sequence. The horizontal axis denotes logarithmically d/Ro, where d is the separation between the two 
stellar centers and Ro the stellar radius at infinite separation. The thick solid and thick dashed lines are reference ones in order 
to check the inclinations of the results easily. 
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Next, we compare the internal velocity field in the co-orbiting frame with that of Taniguchi & Nakamura's analytic 
solution. We focus on the velocity component along the orbital axis (z-axis), because it is three orders of magnitude 
smaller than the x and y-axis components even for the case of closer separation, and is therefore a very valuable 
quantity to check whether the equation of continuity ( |63| ) is well solved or not. In Figs, [l^ - [l5|, we show the velocity 
z-component as a function of the radial distance rn\ from the center of star 1 along three directions: (#m , '/'(l) ) = 
(7r/4,7r/4), (7r/4,7r/2), and (7r/4, 37t/4) and for the orbital separations d = 200 km, 140 km, 100 km, and 70 km. 
It is found that the numerical results agree nicely with those of analytic calculations. Once again, note that the 
discrepancy at small separation comes from the fact that the analytic solution deviates substantially from the exact 
solution. 
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FIG. 13. The z-axis component of the internal velocity field in the co-orbiting frame compared with Taniguchi & Nakamura's 
analytic solution [p3p3|. The horizontal line denotes the radial distance from the center to the surface of star 1, in the direction 
(9n),(fm) = (?r/4, 7r/4). The four panels are snapshots at different separations: 200 km, 140 km, 100 km, and 70 km. Solid 
and dashed lines denote the results of numerical and analytic calculations, respectively. 
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Z-axis component of internal velocity field (9 ^ _ ( 7l /4 j w /2) 
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FIG. 14. Same as Fig. |l3| but for the direction (6(i),<P(i)) = (n/4,n/2). 
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3. Comparison with previous numerical solutions 



As a final test for Newtonian computations, we compare our results with those of Uryu & Eriguchi p3 j for polytropic 
equation of state with 7 = 5/3, 2 and 3, corresponding to polytropic indices n — 1.5, 1 and 0.5 respectively (7 = 
1 + 1/n) . The comparison is presented in Table Q, where the upper lines for each configuration are the results of Uryu 
& Eriguchi and the lower ones are ours. We have chosen the configurations d = 3.6 in Tables 2, 4, and 5 of Uryu & 
Eriguchi j43|. Here, our results are calculated by using N r x Ng x N v = 33 x 25 x 24 spectral coefficients in each of 



the 6 domains and we adopted the same definitions as in Uryu & Eriguchi's article 43 1, namely 

do 
R ' 



d G := (123) 



^(OT< (124) 
J:= (GAPRo) 1 / 2 ' (125) 

^■■=G^JR^ (126) 

where dc is the distance between the two stellar centers of mass and po := M/(4ttRq/3). 

One can see from Table | that our results coincide with those of Uryu & Eriguchi within 0.3 % for physical values 
such as the total energy, the total angular momentum and the orbital angular velocity. Note that the label d of Uryu 
& Eriguchi fli"3"| configurations is the orbital separation between the geometrical centers of two stars normalized by 
the geometrical radius of the star along the x-axis. In our computation, since the geometrical separation is obtained 
after calculation, we cannot fix d initially. Therefore we use the corresponding separation between the centers of mass 
of two stars which Uryu & Eriguchi gave in their paper [|43"f as the orbital separation between the centers of two 
stars d = d /Rp. Although our definition of the center of the star, which is the location of the maximum enthalpy 



(Sect. |VAj), is different from the center of mass, the relative difference between these centers is only about 0.01 % 



around d = 3.6. 



TABLE I. Comparison with the results of Uryu & Eriguchi (1998). 



Separation 






J 


E 


7 = 3 (n = 0.5) 


d G = 3.804 


0.2219 




1.385 


-1.241 


d = 3.804 


0.2211 




1.385 


-1.242 


7 = 2 (n = 1) 


d G = 3.753 


0.2259 




1.371 


-1.133 


d = 3.753 


0.2252 




1.373 


-1.133 


7 = 5/3 (n = 1.5) 


d G = 3.726 


0.2279 




1.364 


-0.9921 


d = 3.726 


0.2274 




1.367 


-0.9911 
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C. Test of the Newtonian limit of relativistic calculations 



We have made many tests of the code in the Newtonian regime up to now, so that we are rather confident in 
the accuracy of the Newtonian part of the code. As a next step, we compare the results of relativistic calculations 
with small compactness (M/R = 7.18 x 1CP 5 ) with those of Newtonian ones. In Figs. [l6| - |l8|, the total energy, 
the total angular momentum and the orbital angular velocity are shown along a sequence. We use N r x Ng x N v = 
25 x 17 x 16 spectral coefficients in each domain. It appears clearly that the results of the small compactness relativistic 
computation coincide with those of the Newtonian computation, as it should be. 
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FIG. 16. Total energy of a relativistic sequence of small compactness (M/R = 7.18 x 10 ) compared with that of that of a 
Newtonian sequence of the same mass. Solid line with filled circles denotes the Newtonian computation and dashed line with 
squares denotes the relativistic one. 
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FIG. 17. Same as Fig. [hi but for the total angular momentum. 
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Orbital ang. velocity compared with Newt, limit of GR calculation 
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FIG. 18. Same as Fig. |l^ but for the orbital angular velocity. 



D. Comparison with previous relativistic numerical solutions 



1. Corotating case 



As a check of for relativistic computations, we compare our results for corotating configurations with those of 
Baumgarte et al. (^6|. We have chosen Tables III and VI of their paper, and compare the results of two different 
separations za — 0.2 and 0.3 in each table. We use N r x Ng x N v = 33 x 25 x 24 spectral coefficients in each domain 
and the criterion 8H = 10~ 7 to end the computation of one configuration. We adopt the same equation of state 
(polytropic with 7 = 2), the same value of the separation rc = d/2 and the same value of the baryon mass Mo. These 
results are shown in Table || where the upper lines for each configuration denote the results of Baumgarte et al., and 
the lower ones correspond to our results. We find a relative discrepancy of 2 % on fj, 4.5 % on g max ; 0.07 % on M, 
0.6 % on J, 4.5 % on 7% and 1.5 % on ¥b- 



TABLE II. Comparison with the results of Baumgarte et al. (1998) 





Mo 


^max 


M 


J 




Q 


TA 


rc 


TB 








TABLE III (M/R = 


0.05) 










0.20 


0.0595 


0.0284 


0.057815 


0.01109 




0.048 


0.591 


1.791 


2.959 






0.0280 


0.057816 


0.01113 




0.048 


0.582 




2.923 


0.30 




0.0288 


0.057836 


0.01155 




0.037 


0.975 


2.118 


3.251 






0.0285 


0.057836 


0.01161 




0.038 


0.968 




3.217 








TABLE VI {M/R = 


0.15) 










0.20 


0.1534 


0.1303 


0.140859 


0.04174 




0.116 


0.413 


1.244 


2.067 






0.1242 


0.140774 


0.04194 




0.117 


0.395 




2.037 


0.30 




0.1341 


0.140971 


0.04268 




0.092 


0.682 


1.477 


2.273 






0.1286 


0.140874 


0.04294 




0.093 


0.668 




2.249 



2. Irrotational case 



For irrotational relativistic configurations, a detailed comparison with Uryu & Eriguchi results [j23|,|24| is underway 
p5[ . For the purpose of the present article, we have compared only the cusp point configuration of a M/R = 0.14 
7 = 2 polytropic sequence as given in the last but one line of Table IV of Ref. |^3| and the last line of our Table IV 
below. The agreement is very satisfactory: the relative discr epancy is below 0.1% on M, 0.4% on QMb, 2.1% on 
J/M 2 and 2.1% on the semi-diameter ap of the stars [Eq. |129|) below]. 
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VI. RESULTS FOR 7 = 2 POLYTROPES 



The tests of the code being successfully p asse d, we employ the code to compute an irrotational relativistic sequence 
based on the polytropic equation of state (105)-(107) with 7 = 2. We are using n — 0.0332 ( o nu 1 c c 2 , and consider the 
compactness parameter M/R = 0.14 at infinite separation. This results in the baryon mass Mb = 1.625M©. 

The computational parameters are as follows: six domains are used, such that (using the notations of Sect. IV A| ): 



iV/i\ = N{2) = 3, Mn\ — M(n\ — 1, with the same number of coefficients in each domain: N r x NgX N v = 33 x 21 x 20. 
The criterion to end the computation of one configuration is set to 5H = 10 -7 . A special treatment has been performed 



for the closest configuration, because of the existence of a cusp on the stellar surface (Sect. IV E): N r x Ng x N v = 
25 x 17 x 16 coefficients have been used along with the enthalpy gradient threshold Xb = 0.55, resulting in a frozen 
mapping. 

In Fi gs. O and ||(], the half of the ADM mass and the total angular momentum of the binary system, as defined 
in Sect. |lIID[ arc shown along the evolutionary sequence. This sequence ends at around d = 37.5 km (/ = 380 Hz) 
where a cusp appears on the surface of the stars. 

One can see from these figures that there is no turning point for the 7 = 2 case. This result agrees with that of 
Uryu & Eriguchi pi. 
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FIG. 19. Half of the ADM mass of the binary system as a function of the coordinate separation for an evolutionary sequence 
of relativistic irrotational stars. 
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FIG. 20. Same as Fig. [Isj but for the total angular momentum. 
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An important result of this computation has already been presented in Ref. plfl , namely the central energy density 
remains rather constant (with a slight increase below 0.01 %) and finally decreases (see Fig. ^|). As discussed in the 
Introduction, this result makes the collapse of the individual neutron stars to black hole very unlikely prior to the 
merger. 



We summarize the results about the sequence in Table III, listing the ADM mass, the total angular momentum, the 



orbital angular velocity, the axis ratios, and the relative change in central energy density along the quasiequilibrium 
sequence. Since this is the first table presented for a sequence of relativistic irrotational binary neutron stars, we give 
a rather large number of digits in order to compare with the results of other works from now on. Note that we are 
using the following values of the fundamental constants: G — 6.6726 x 10 -11 m 3 kg _1 s -2 , c = 2.99792458 x 10 8 ms _1 
and M Q = 1.989 x 10 30 kg. 

TABLE III. Half of ADM mass M, total angular momentum J, orbital angular velocity f2, axis ratios, and relative change 
in central energy density along a Mb = 1.625 Mq quasiequilibrium sequence constructed upon a 7 = 2 polytropic EOS. <Xi, 
<i2, and 03 denote the coordinate lengths parallel to the semi-major axes x, y, and z, respectively. ai, op p is the length in the 
direction opposite to the companion star. 



[km] 


0.5 x M [M Q ] 


J [GMl/c] 


[rad/s] 


n/(2vr) [Hz] 


02/01 


03/01 


Ol,opp/Ol 


(6 C 6c, 00 ) /Cc,oo 


100 


1.50545 


11.8370 


597.24 


95.054 


0.99100 


0.99367 


0.99319 


4.0606e-05 


90 


1.50457 


11.3403 


695.15 


110.64 


0.98839 


0.99139 


0.99234 


6.0695e-05 


80 


1.50351 


10.8243 


823.17 


131.01 


0.98445 


0.98788 


0.99106 


6.9369e-05 


70 


1.50223 


10.2880 


996.14 


158.54 


0.97811 


0.98210 


0.98903 


8.4666e-05 


60 


1.50065 


9.73115 


1239.9 


197.34 


0.96687 


0.97171 


0.98466 


3.2735e-05 


50 


1.49870 


9.15576 


1603.5 


255.20 


0.94402 


0.95037 


0.97369 


-5.9816e-04 


45 


1.49758 


8.86296 


1858.8 


295.84 


0.92226 


0.92999 


0.96263 


-2.0684e-03 


42 


1.49679 


8.68172 


2041.4 


324.89 


0.90315 


0.91100 


0.95408 


-4.2246e-03 


41 


1.49655 


8.62425 


2111.0 


335.98 


0.89206 


0.89940 


0.95076 


-5.4239e-03 


37.5 


1.49572 


8.43623 


2389.7 


380.33 


0.81445 


0.82752 


0.91949 


-1.2238e-02 
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For comparison purposes with other works, we also give Table [fv| in which the physical quantities are normalized 
by using the equation of state constants k and 7 to set a length scale i? po i y according to 

(127) 



#poly := k 2 ^" 1 ' ■ 

We can therefore define the same dimensionless quantities as Baumgarte et al. [^6| and Uryu & Eriguchi [£3j : Mb = 
AW-Rpoiy, M = M/R po \y, J = J/i?^ oly , = R po \yft, d= d/Rpoiy, d G = d G /R po i yi a = a /R po i y , where d G is the 
distance between the "centers of mass" of each stars as defined by Eq. (107) of Uryu & Eriguchi p3): 



A 3 T D nXd 3 x - 



star 1 



NL 



(2) 



A 3 T n nXd 3 x 



star 2 



(128) 



(129) 



and ao is half of the coordinate length of a star along the X axis: 

^0 — ~ I ^max ^min | • 

This latter quantity is denoted i?o by Uryu & Eriguchi [^3| . 

TABLE IV. Dimensionless ADM mass M, total angular momentum J, orbital angular velocity Q , an d half of the coordinate 
length along the X axis ao along the Mb = 0.146202 quasiequilibrium sequence presented in Table 



III 



d 


d,G 


0.5 x M 


J 


J/M 2 


Q. 




ao 


6.0927 


6.0924 


0.135446 


9.58174e-02 


1.30573 


3.2698e-02 


4.7805e-03 


0.81089 


5.4835 


5.4830 


0.135367 


9.17965e-02 


1.25239 


3.8058e-02 


5.5642e-03 


0.80934 


4.8742 


4.8736 


0.135272 


8.76195e-02 


1.19708 


4.5067e-02 


6.5889e-03 


0.80770 


4.2649 


4.2642 


0.135156 


8.32782e-02 


1.13973 


5.4536e-02 


7.9733e-03 


0.80623 


3.6556 


3.6546 


0.135014 


7.87708e-02 


1.08031 


6.7883e-02 


9.9246e-03 


0.80542 


3.0464 


3.0447 


0.134839 


7.41131e-02 


1.01907 


8.7788e-02 


1.2835e-02 


0.80796 


2.7417 


2.7396 


0.134738 


7.17429e-02 


0.98796 


1.0177e-01 


1.4879e-02 


0.81489 


2.5589 


2.5565 


0.134667 


7.02758e-02 


0.96878 


1.1176e-01 


1.6339e-02 


0.82583 


2.4980 


2.4954 


0.134645 


6.98107e-02 


0.96268 


1.1557e-01 


1.6897e-02 


0.83271 


2.2848 


2.2810 


0.134571 


6.82887e-02 


0.94273 


1.3083e-01 


1.9127e-02 


0.88016 
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Finally, let us show some figures about metric, density, and internal velocit y fiel ds. The lapse function N is 
represented in Fig. ^l|. The coordinate system is the (X, Y, Z) one defined in Sect. IV A and the coordinate separation 
is 41 km (last but one line in Tables III and IV). At this separation, the central value of N of each star is 0.6416. 

The shift vector N of non-rotating coordinates [defined by Eq. (B7f)] is shown in Fig. ^2|. The plot is a cross section 
of the orbital plane. 

The K xx , K , and K YY components of the extrinsic curvature tensor of the hypersurfaces t = const are shown 
in Fig. ^3|. The values in the figures are multiplied by the square conformal factor A 2 , and the plots are cross section 
of the orbital plane. 



ln(N) (x=xl) 



ln(N) (y=0) 



ln(N) (z=0) 




y [km] x [km] x [km] 

FIG. 21. Isocontour of the gravitational potential v (logarithm of the lapse function N) when the coordinate separation is 
41 km. The plots are cross sections of the X = —20.5 km, Y = 0, and Z — planes. The thick solid lines denote the surfaces 
of the stars. 



Shift vector (z = 0) 




x [km] 

FIG. 22. Shift vector N of non-rotating coordinates in the orbital plane when the coordinate separation is 41 km. The 
thick solid lines denote the surfaces of the stars. 
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A 2 K D (z = 0) 



A 3 K xy (z = 0) 



A 2 K yy (z=0) 




x [km] 



x [km] 
YY 



FIG. 23. Isocontours of the components K x A , K XY , and K Y (multiplied by A 2 ) 
the coordinate separation is 41 km. The plots are cross sections of the orbital plane (Z ■ 
positive (resp. negative) values. The thick solid lines mark the surfaces of the stars. 



x [km] 

of the extrinsic curvature tensor when 
= 0). Solid (resp. dashed) lines denotes 



The baryon density field in the fluid frame is shown in Fig. The plots are cross sections of Z — and Y = 
planes. 

Finally, we show in Fig. |25| the internal velocity field in the co-orbiting frame, or more precisely the orthogonal 
projection in the S t hypersurface of the vector field V given by Eq. ([30]). Note that this vector field is tangent to the 
surface of the stars, as it should be. 

Baryon density (z = 0) Baryon density (y=0) 



oo 



. L i i i i I i i i i J 

-50 50 '-50 50 

x [km] x [km] 

FIG. 24. Isocontours of the fluid proper baryon density n when the coordinate separation is 41 km. The plots are cross 
sections of Z = and Y — planes. The thick solid lines denote the surfaces of the stars. 
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Velocity w.r.t corotating frame (z = 0) 
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FIG. 25. Internal velocity field with respect to the co-orbiting frame in the orbital plane when the coordinate separation is 
41 km. The thick solid lines denote the surfaces of the stars. 

VII. DISCUSSION 
A. Comparison with other numerical methods 

The numerical method presented in this article is the only method for computing relativistic binaries in which the 
computational domains extends to infinity, thereby enabling to impose exact boundary conditions on the gravitational 
field equations. All the other methods p6| , p2|]25| ] employ finite computational boxes. Our experience from calculations 
of single rotating neutron stars show that the finite size of the computational domain can result in some loss of 
accuracy (see Ref. Q for a discussion of this point). 

Besides, thanks to the splitting of the metric quantities in a part described on the domains centered on star 1 and 
a part described on the domains centered on star 2, we can describe without any loss of accuracy very distant stars. 
In fact we can recover the spherical limit when the stars have very large separations, contrary to all other numerical 
methods which are losing resolution when the stars are put farther apart (see for instance the discussion in Sect. VA 
of Ref. H). 

We are using surface-fitted spherical coordinates, which by construction are well adapted to describe the stellar 
fluid interiors. As it can be seen on Figs. ^lf]2li| , these coordinate systems, which are centered on one of the two stars, 
are also well adapted to the description of the metric quantities, because these latter are maximum at the location 
of the stars. Baumgarte et al. |^6| and Marronetti et al. |22| use instead Cartesian coordinates in a single domain 
("box"). Closer to our approach, Uryu & Eriguchi |23| developed a multi-domain method with surface-fitted spherical 
coordinates, which enable them to treat precisely the fluid interiors of the stars. However, for the gravitational field 
they use a single spherical coordinate system which is centered at the system center of mass. 

As far as irrotational binaries are concerned, we payed a special attention to the resolution of the equation for the 
velocity potential First we solve numerically only for a small part 'Fo of 'F, thereby reducing the numerical error. 
Second we let appear in the equation for 'Fo a partial differential operator which is invertible and give as a unique 
solution that with the correct behavior at the stellar surface (velocity field tangent to the surface in the co-orbiting 
frame). The equation for 'F is instead solved as a Poisson equation A'F = source with a boundary condition at the 



stellar surface by Uryu & Eriguchi |23J] and Marronetti et al. [22|31J]. Note that Marronetti et al. performs only an 
approximate treatment of the boundary condition, which amounts to considering that the surface of the star is an 
exact sphere. This is of course not valid for close configurations. On the contrary, thanks to the introduction of 
surface-fitted coordinates, Uryu & Eriguchi |^3| have been able to treat the boundary condition exactly. 

B. Tests passed by the code 

We have performed extensive tests of the numerical code. In particular, we have shown that, in the Newtonian 
limit, our numerical results coincides with the semi- analytical solutions recently obtained by Taniguchi & Nakamura 
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p3| , p4| for compressible polytropic stars. The only discrepancies appeared to be due to missing higher order terms in 
Taniguchi & Nakamura's solutions and not to some inaccuracy of the numerical code. 

Regarding relativistic configurations, no analytical solutions was available to compare with. In this case, we checked 
only by comparing with previous numerical solutions, namely that of Baumgarte et al. |B6[ for synchronized binaries 
and Uryu & Eriguchi |23| for irrotational binaries. The agreement is of the order of 1%. For the astrophysically 
relevant case of irrotational relativistic binaries, a detailed comparison with Uryu & Eriguchi code [^3| is underway 
& 

C. Future prospects 

We are currently using the method described in this article to compute models of close binary neutron stars with 
various equations of states: polytropic EOS with various polytropic indices, dense matter EOS resulting from recent 
nuclear physics calculations. In particular, we are studying how parameters like the frequency location of the innermost 
stable orbit (if any) depends on the equation of state, in order to help in the interpretation of gravitational wave 
signals from coalescing neutron star binaries. The results of these studies will be published elsewhere |67]]. 
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APPENDIX A: LINK WITH TEUKOLSKY'S AND SHIBATA'S FORMULATIONS 

The first integral of motion for quasiequilibrium irrotational binaries derived by Teukolsky |]l9| is (cf. his Eq. (57), 
re-written within our notations^) 

7V(/i 2 + D1'-D*) 1/2 + B-D* = C , (Al) 

where C is some constant. At first glance, this might look quite different from our first integral of motion (|34|). 
However, if we substitute Eq. (|2^) for T in the exponential form of our first integral (i.e. Eq. (|33|)) we get 

hNT n (l - U • U ) = const. (A2) 

By means of Eqs. (|3l]) and (Q), this equations becomes 

WVT„ + B = const. (A3) 



Finally, if we substitute Eq. (J33) for T n is this relation, we recover Eq. (Al). In particular this shows that the constant 
in the right-hand side of Eq. (|33|) is nothing but the constant denoted C by Teukolsky [ pi . 

The first integral of motion for quasiequilibrium irrotational binaries derived by Shibata |2Ct| is (cf. his Eq. (2.18), 
re-written within our notations^) 

^ + S • D* = const, (A4) 
A 

where A and S are defined by the following decomposition of the fluid 4-velocity in a part along the Killing vector 1 
and a part parallel to the hypersurface E t : 

u = A(l + S) with n • S . (A5) 



Our shift vector B is the negative of Teukolsky 's B. 
9 Our A and S are denoted respectively by u° and V by Shibata 
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Now, substituting Eq. (||) for 1 in this relation and using Eq. (|27j), we get 



S = ir n u + B= -5-D* + B 

A Ah 



(A6) 



where the second equality follows from Eq. (pl|). Inserting Eq. (A6) into (A5) and using the normalization relation 
u • u = —1 results in the following expression for A: 



1 

hN 



A = (h 2 + • D* 



1/2 



(A7) 



Finally substituting Eq. (A7) for A a nd E q. ( |A6| ) for S into Shibata's first integral of motion (A.4) result in Teukolsky's 
form of the integral of motion (Eq. ( |Al|) above), which shows the equivalence of the various formulations. 



APPENDIX B: NUMERICAL METHOD TO SOLVE THE ELLIPTIC EQUATION FOR THE VELOCITY 

POTENTIAL 



Equation ( p3|) for the part \I/o of the velocity potential can be written 
with 

a := (H , 



(Bl) 
(B2) 



and 



b l := (1 - CH)VH + CHVp 



n := [W - n7,)V,// Jl [ WfoiiH-fl + ^ViTn 



(B3) 



(B4) 



Equation (Bl) is not merely a Poisson type equation for ty because the coefficient a vanishes at the surface of the 
star. It therefore deserves a special treatment. In the works of Marronetti et al. p2|,|3l]l and Uryu & Eriguchi p3[ , 
Eq. (Bl) is recast as a Poisson equationp°| A^>o = source, dividing both sides of Eq. (|l3l|) by a. In order that the 
right-hand side be regular, one must then impose as a boundary condition fe l Vi*o — cr = at the surface of the star. 
We choose here a different approach, considering that the operator in Eq. (Bl) is not the Laplacian but instead the 
operator 



LV :=a(l-£ 2 )A ? *o + /3£ 



(B5) 



where a and (3 are two constants, £ 6 [0, 1] is the computational radial coordinate introduced in the mapping ([74]), 
and is an operator which, expressed in terms of the computational coordinates (£,#',<//) has the same structure 
than the Laplacian operator: 



A^o := 



1 9 



9* 



9 ( . ,,9*0 
sin 



£ 2 <9£ V S di ) ' £ 2 sin 0' dO' V""" d9' J ' sin 2 6' dip' 2 



9 2 * 



(B6) 



Here we assume that there is only one domain covering the star, i.e. that Mm\ — M{ 2 ) = 1, so that the surface of the 
star is given by £ = 1. Equation (Bl) is then re- written as 



L*o = a + a(l- C 2 )A c *o - aAfo + P£ 



9*o l4 = 



9£ 



(B7) 



3 These authors are using *!/ and not as the unknown function, but this has no consequence on the following discussion. 
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The basic idea is to solve this equation by iterations, considering the whole right-hand as a source term and using the 
previous step value of Wo in it- One must also choose the constants a and f3 so that the term a(l — ^ 2 )A^q [resp. 
P^d^o/dS] is as close as possible to aMo [resp. b l VjWo]. We opt for the following choices: 



and 



where b r is the r-component of the vector b l . In solving Eq. 



(B8) 



(B9) 



by iterations, we introduce the following relaxation 



\*Z + (i - \)* J 



j -i 



(BIO) 



where J (resp. J — 1) labels the current step (resp. previous step) and A is the relaxation factor, typically chosen to 
be 0.5. 

At each step of the iteration, Eq. (B7) is solved by the following spectral method. First an expansion in spherical 
harmonics Y" l (9',ip') is performed, so that Eq. ( p37| ) becomes equivalent to a set of ordinary differential equations 
[one equation for each couple (£, m)]: 



(Bll) 



where /(£) and s(£) are the (£,m) coefhcient of Wo ancl of the whole right-hand side of (B7) respectively and Lg is 
the following differential operator: 



L e f:=a(l-e) 



d 2 f 2df £{i+l) 



(B12) 



Since the source s(£) vanishes for t = Q, we treat only the case £ > 0. Regularity properties at the origin (£ = 0) 
imply that /(£) and s(£) should be expandable in even (resp. odd) Chebyshev polynomials T n (£) for £ even (resp. 
odd). Due to the division by £ and £ 2 , the differential operator Lg is singular on Chebyshev polynomials at £ = 0, 
except for t = 1. Therefore, instead of Chebyshev polynomials, we use the following polynomials P n (£) as a expansion 

basis for / (N is the total number of coefficients in the Chebyshev expansions, denoted 7vj Q> (0) in Sect. [IV B|) 
. for * even : P n (£) := T 2n (0 + T 2n+2 (0 = 2£T 2n+1 (£), < n < N - 2; 

• for i = 1 : P„(0 := T 2n+ i(0, < n < JV - 1; 

• for ^ odd and £ > 1: P n (£) := (2n + 3)T 2 „+i(£) + (2n + l)T 2 „ +3 (0, < n < - 2. 



The operator on the left-hand side of Eq. (Bll) is regular for each of the polynomials P n (£) (such a basis is called a 
Galerkin basis). 
We thus consider the differential operator Lg acting 

• for £ even : from the N — 1 dimensional vectorial space span by the polynomials P n (£) (0 < n < N — 2) to the 
N — 1 dimensional vectorial space span by the polynomials T 2n (C) (0 < n < — 2); 

• for £ = 1 : from the N dimensional vectorial space span by the polynomials P n (£) = P2n+i(£) (0 < n < N — 1) 
to itself; 

• for £ odd and £ > 1: from the A — 1 dimensional vectorial space span by the polynomials P n (0 (0 < n < A — 2) 
to the N — 1 dimensional vectorial space span by the polynomials T 2n +i(£) (0 < n < N — 2). 

The operator Lg is then one-to-one (isomorphism) between these vectorial spaces. T his m eans that the only homoge- 
neous solution is zero. Otherwise stated, for each I there is a unique solution to Eq. (Bll) in the vectorial space spans 
by the P n (£)'s. To find this solution, we transform the matrix Aij of Lg in the above bases into a banded matrix by 
means of the following linear combinations: 
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• for £ even 



— 



i + V 



[A i:j - A {i+1)j ] for < i < N - 3 



A^ = Aij - A (l+2)J for < i < N - 5 



(B13) 
(B14) 



• for £ odd 



Ai 



.4 



(*+2) 



j] for < i < N - 3 ; 



A^ = A i3 ~ A {l+2)j for < i < N - 5 



(B15) 
(B16) 



Since the resulting matrix Aij has at most 5 bands, the linear system is easily and CPU-efficiently solved to get the 
coefficients of the solution / in the basis of the polynomials P n (£)- A simple combination is then performed on these 
coefficients to get the coefficients on the usual Chebyshev bases. 

A very discriminating test of this numerical technique, namely the evaluation of the tiny z-component of the velocity 



field resulting from VP, is presented in Figs. pT3|-|T5| (Sect. VB2) 
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